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Abstract. The increased availability of radio occultation 
(RO) data offers the ability to detect and study turbulence 
in the Earth’s atmosphere. An analysis of how RO data can 
be used to determine the strength and location of turbulent 
regions is presented. This includes the derivation of a model 
for the power spectrum of the log-amplitude and phase fluc- 
tuations of the permittivity (or index of refraction) field. The 
bulk of the paper is then concerned with the estimation of the 
model parameters. Parameter estimators are introduced and 
some of their statistical properties are studied. These estima- 
tors are then applied to simulated log-amplitude RO signals. 
This includes the analysis of global statistics derived from a 
large number of realizations, as well as case studies that illus- 
trate various specific aspects of the problem. Improvements 
to the basic estimation methods are discussed, and their ben- 
eficial properties are illustrated. The estimation techniques 
are then applied to real occultation data. Only two cases 
are presented, but they illustrate some of the salient features 
inherent in real data. 


1 Introduction 

There is a long and distinguished history in the study of elec- 
tromagnetic (EM) wave propagation through random media 
(cf. Tatarskii, 1971; Yeh and Liu, 1982; Ishimaru, 1997). 
These works have provided a firm theoretical foundation for 
estimating statistical properties of the neutral atmosphere and 
ionosphere via the statistical properties of the received EM 
signals. That is, the characteristics of the turbulent atmo- 
sphere can be deduced from, for example, correlation and/or 


spectral analysis of the phase and/or amplitude of the re- 
ceived signal. In the past, the bulk of the experimental analy- 
sis in this area has been performed with ground-based trans- 
mitters and receivers (e.g., radars and lidars), as well as with 
ground-based receivers and space-based transmitters. With 
the advent of Global Navigation Satellite System (GNSS) 
constellations, a new avenue has become available to investi- 
gate the turbulence properties of the earth’s atmosphere. The 
deployment of an ever-increasing number of Low Earth Or- 
biting (LEO) satellites - with high quality, high-sampling 
rate receivers - provides a very valuable new source of 
turbulence measurements: GNSS-LEO occultations. 

Previous efforts to study turbulence with RO signals have 
been focused on the study of the middle to upper strato- 
sphere and the ionosphere (e.g., Gurvich and Brekhovskikh, 
2001; Gurvich and Kan, 2003; Gurvich and Chunchuzov, 
2005). On the other hand, very few works have addressed 
the characterization of turbulence in the upper troposphere 
and lower stratosphere. Numerous other practical applica- 
tions can benefit by having turbulence measurements (and re- 
sultant climatologies) from GNSS-aircraft occultations. For 
example: space weather (Hajj et al., 2000), transionospheric 
communication links (Secan et al., 1997), aviation safety 
and navigation systems (Comman et al., 2004; Conker et 
al., 2003), the accuracy of ground-based GNSS receivers 
(GangulvBHiet al.?, 2004), as well as the effect of turbu- 
lence on measuring atmospheric state variables from GNSS- 
LEO occultations (Gorbunov and Kirchengast, 2005). 

In this paper, a detailed analysis of the problem is pre- 
sented, including the derivation of a model for the power 
spectrum of the log-amplitude and phase fluctuations of the 
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permittivity (or index of refraction) field. The bulk of the pa- 
per is then concerned with the estimation of model parame- 
ters, such as the intensity and location of the turbulence along 
the line of sight. A thorough study of the properties of the 
parameter estimators is presented using simulated data. This 
includes global statistics over a large number of realizations, 
as well as case studies that illustrate various specific aspects 
of the problem. The application of the methodologies is then 
made with real occultation data. As the focus in this paper 
was the development of methods, the purpose of the real data 
analysis is to show how the theoretical and practical concepts 
carry over into the real data. Hence, only two cases are pre- 
sented, but they illustrate some of the salient features inher- 
ent in real data. A more thorough analysis of real data is left 
to follow-on efforts. 

2 Wave propagation theory 
2.1 Assumptions 

In order to make the problem tractable, a number of simpli- 
fying assumptions will be made. Most of these are routinely 
used in the literature for optical wavelengths (cf. Tatarskii, 
1971; Ishimaru, 1997). Very little work has been done to 
justify many of these assumptions for the decimeter wave- 
lengths in the GNSS problem; however, a few excellent pa- 
pers exist, including Strohbehn and Clifford (1967), Clifford 
and Strohbehn (1970), Clifford and Lataitis (1985), and Lee 
and Harp (1969). We will list all of the assumptions here, 
but their use will be introduced below in the context of the 
development. 

1. Polarization effects are ignored. This allows for the 
wave equation to be written in terms of a single 
scalar component of the electric field vector. This 
is a reasonable assumption for the decimeter wave- 
lengths consider for the GNSS problem (cf. Clifford and 
Strohbehn, 1970). 

2. Weak and single scattering is assumed. This is rea- 
sonable for the GNSS transmitter wavelengths and 
typical path lengths through turbulence (Clifford and 
Strohbehn, 1970). 

3. The assumption that the transmitter wavelength, A, is 
much smaller than the smallest fluctuations in the per- 
mittivity field, Iq, (A <£ k), allows for a number of sim- 
plifications in the development (including the first two 
items above). While this is not necessarily true for 
decimeter wavelengths, Clifford and Strohbehn (1970) 
show that the practical effects of the assumption A <<C lo , 
can be extended into regime, A > lo, but not A /o. 

4. Taylor’s Hypothesis is valid. That is, for short time pe- 
riods advection dominates over internal fluid stresses 
(Monin and Yaglom, 1972). This is important for a 


medium that is moving, which in general is the case for 
the atmosphere. The practical effect of this assumption 
is that we can assume that the permittivity for a small 
scattering region is constant as it moves during a short 
time interval. 

5. The Markov Approximation is valid (Tatarskii, 1971, 
Chapter 5, Section B). This means that there are mini- 
mal effects on the received signal due to correlation in 
the permittivity field along the line-of-sight. Therefore, 
one can assume that this correlation is represented by a 
delta-function. 

6. The transmitter is a point source (i.e., we do not con- 
sider the radiation pattern), and the receiver acts as point 
aperture along the line-of-sight (i.e., we do not consider 
the gain pattern). Due to the typical large distances be- 
tween the transmitter and the scattering medium, the 
first assumption is quite reasonable. If the distance be- 
tween the scattering volume and the receiver is also 
large, and a high-gain antenna is used, the second as- 
sumption is also a reasonable one. For an omnidirec- 
tional receiver antenna and/or the receiver very close to 
the scattering volume (e.g., for an airborne receiver in 
the atmosphere), then these effects do need to be ac- 
commodated (cf. Lee and Harp, 1969; Tatarskii, 1971, 
Section 53). 

7. Atmospheric attenuation due to absorption is assumed 
to be smooth and varying slowly during an occulta- 
tion, and so that its effects can be accommodated by 
trend removal of the amplitude or phase signal in the 
time domain. This also means that we can consider the 
permittivity field to be real-valued. 

8. The turbulence patch has compact support. This is re- 
quired to ensure the existence of various integrals that 
arise in the calculations. 

9. The effects due to an inhomogeneous background (e.g., 
refraction) are minimal. This means that we can assume 
straight-line propagation. 

10. The medium is purely random. There are no 
deterministic structures, such as layered phenomena. 

Most of the assumptions presented above are reasonable for a 
first-order analysis of the turbulence problem. Nevertheless, 
there are two aspects of the problem that are not addressed 
in this initial study, the effects due to a background or turbu- 
lent permittivity field characterized by inhomogeneity and/or 
anisotropy. These are obviously real-world situations, but as 
we are trying to keep to a simple analysis at this point, they 
are not included in what follows. 
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2.2 The first Rytov approximation 

If polarization effects can be ignored, the vector wave equa- 
tion for monochromatic signals, derived from Maxwell’s 
equations can be written in the scalar form, 

(V 2 + k 2 s(r))U (r) = 0 (1) 

where U is a component of the electric field (a transverse 
component for our case), k is the (constant) wavenumber for 
the transmitted frequency, k = Ircf/c = 2n/k, and e is the 
permittivity of the medium. The index of refraction, n , is 
given by £ = n 2 , and free space is characterized by £ = 1. 
Since we have assumed weak scattering, a perturbation ap- 
proached is used. Re-write Eq. (37) as 

(V 2 +k 2 )U (r) = —k 2 e'(r)U(r) (2) 

where e = {s)+e', and s' <3C ( e ). It is assumed that (e) = 1, 
its ffee-space value. The angled brackets refer to an ensem- 
ble average. Standard approaches to solving this equation 
are provided in the literature (e.g., Tatarskii, 1971; Ishimaru, 
1997) - one of which is the first-order Rytov approach. We 
return to the Helmholtz wave equation and let 

U(r) = e ;* (r) (3) 

where the complex phase, Vc is given by 

V Hr) = X (r) + iS(r) 



Fig. 1. Scattering geometry at time /] . 


we will want the origin of the coordinate system to be at the 
(local) center of the earth. In the arbitrary coordinate sys- 
tem, the receiver is at the point /*r, the transmitter is at rj, 
and a given scatterer is at r $ (see Fig. 1) (we are consider- 
ing a continuum field of permittivity, so the term scatterer is 
used loosely). We assume that multiple scattering does not 
occur, and that the transmitter acts as a point source, (i.e., an 
outgoing spherical wave), which leads to 


and x is the log-amplitude function, and S is the phase 
function. Plugging Eq. (39) into Eq. (37) and solving 
for C/oV'i leads to the first-order term (cf. Tatarskii, 1971; 
Ishimaru, 1997), 


* ,{r) =v^f 

Where 


Go(r-r')£'(r')Uo(r')dr' 


Go(r-r'): 


II II 

4 n II r-7\\ 


( 4 ) 


( 5 ) 


HHL/ = differential operator? is the free-space Green’s 
function for the problem (the double vertical lines indicate 
the vector magnitude). Equation (40) is known as the first 
Rytov approximation. 


2.3 Correlation function and frequency spectrum 


We depart from the standard development as presented by, 
for example, Tatarskii (1971) or Ishimaru (1997), going for- 
ward. We consider the situation where the transmitter, re- 
ceiver, and media are all moving relative to each other, as 
well as with respect to a fixed coordinate system. Typically, 
the origin of the coordinate system is set at the transmitter 
or at a point within the scattering volume. However, in what 
follows it is assumed that the origin of the coordinate sys- 
tem is at an arbitrary point. For the GNSS-LEO geometry, 


Vm (*"r) = IkR — r T ||c _i * l|rR — rT " 

An 

/ ^/<:||r R -r s || g i/:||rs-r T || 

z | 7 T -s'(r s )drs (6) 

IkR-rsll Iks -n- 1| 

In order to ensure the existence of the integral, we assume 
that the permittivity fluctuations have compact support, and 
s'(rs) is non-zero (or has a minimal effect on Vm(7*r)) in a 
region close to the line of sight, r R — rj. If £'(r$) is a random 
field, we assume that Eq. (42) exists in the sense that it is a 
member of an ensemble of random permittivity fluctuations. 
We will assume straight-line propagation. Return to Eq. (42). 
If the scattering angles are small - a reasonable assumption if 
the ratio \/Iq is small - we can make the so-called parabolic 
approximation in the integral equation (Tatarskii, 1971, Sec- 
tion 45) (/o is a length scale on the order of the smallest fluc- 
tuations in the medium). In order to use this approximation, 
we will consider basis vectors that are parallel or perpendic- 
ular to the vector r r — rj. We will denote the component of 
vectors parallel to tr — rp, “jc” and the perpendicular parts, 
p (i.e., the projection of a given vector onto the plane per- 
pendicular to r r — ry). Note that p R — p T = 0. Consider 

the quantity, IkR-^sll = [Ur-^s) 2 + ||/> r - / o s ||] I/ '. We 
will assume that xj < Jts < x R . If jcr — jts, the distance along 
the x axis from the scatterer to the receiver, is much greater 
than the transverse distances (valid for small-angle scatter- 
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ing), then ||tr — rsll, can be approximated by 

,,2-11/2 


lkR-rsll = (JCR-*s) 


s *R-*S + 


1 + 


IPR-PS 


(XR~X S ) 

l<°R~/ ) sll 

2(xr-x s ) 


(7) 


This is the parabolic approximation. A similar expan- 
sion can be made for ||/*r — rj\\, which is simply *r — xj. 
Since k ||tr — r s|| = 2n || r r — rsll A 1, we will use the 
full expression in Eq. (43) for the exponential terms in 
Eq. (42); whereas, we only need to keep the leading term for 
the non-exponential terms. Plugging all these factors back 
into Eq. (42), we have for the non-exponential terms 

IkR— r T || _ *R--*T , g . 

IkR -rs II Iks - '•-rll (-XR -xs) ( XS - XT) 


which we will denote as b(x&. jcs,*t)- It can be shown that 
the exponential term from Eq. (42) is given by 


IkR -rsll + Iks -'•tII - IkR — /•tII 

b(XR,Xs,X T ) |i 2 

% ~ |/>s-pt| 


(9) 


Plugging Eqs. (44) and (45) into Eq. (42) gives the desired 
result, 


VMk r) 
exp 


-I 

4 IT J 


b(XR.XS,XT) 


b(x K .x s.xt),| |,2l ,, ,, 

ik ||/ 0 s _ /°t|| £ ( r s)dr$ 


( 10 ) 


If we move the origin to rj, this is equivalent to setting xj = 
0 and p T = 0, and since we have assumed that p R — p T = 
0, this means that p R = 0. Equation. (46) then describes 
the field, 


tl(XR'PR 


k 2 f 


*S (XR-X$) 


exp 


ik- 


xr 

2xs {xr-xs) 


ll^sll 


e'(r s )drs 


( 11 ) 


which is equivalent to Eq. (49) — (5) in Tatarskii (1971). 
Therefore, we see that Eq. (46) has the correct limiting form 
when the origin is placed at the transmitter. 

Next, consider the transmitter, receiver and atmosphere 
moving relative to each other - and all with respect to a 
fixed origin, e.g., the center of the Earth. Since the longi- 
tudinal motions do not change the value of the correlation 
of the field in an appreciable manner (Klyatskin, 1970), we 
assume that the ’'-positions are not changing, and so the 
velocity vectors described below are the projection of the 
three-dimensional velocity vectors onto the plane perpen- 
dicular to /*r — rj. We consider a short enough time in- 
terval, t = ?2 — 1\, such that Pi{ti) = Piit\) + Vi t, (where 
“z” stands for R, T, or S). For this work we will assume 



Fig. 2. Scattering geometry at times t\ and ti- 


that the velocities are constant over all (short) time intervals 
(i.e., Vi( r) = Vi); that they are the same for all the scatter- 
ed (i.e., Vs, = Vs,); and that they are deterministic quan- 
tities (i.e., (Vi) = Vi). The same procedure used to derive 
Eq. (46) can be applied to calculating ^\(r R (t 2 )), except 
that we now use basis vectors parallel and perpendicular to 
rR — ri + (VR— Vj)r (see Fig. 2). Unfortunately, this will 
complicate matters later on when we calculate the correla- 
tion functions; which in turn contain products of \j/\ (rR(fi)) 
and yjr\ (r R (t 2 )). That is, we will have to transform the com- 
ponents of vectors in one basis into the other. However, if 
the two vectors r R — rj and /*r — tt + (Vr — ViTrare close 
to parallel we can then choose the basis vectors perpendicu- 
lar to them to be identical. This assumption is valid to first 
order if ||V R r|| « ||r R -r T || and || V T r II « lk R -r T ||. In 
this case, the components of a vector in each basis system 
will be the same - up to the level of approximation used. 
Given the distances involved in typical GNSS-LEO occul- 
tations, these conditions will be met; hence, we can safely 
assume that the same basis can be used for calculations in- 
volving the components of both r R — tt-KVr — VtK and 
r K -r T . 

Next, we calculate \f/\ (rR^)). Consistent with the 
approximations just discussed, we assume that 

*R-*s» ||/0 r -/0 S + (Vr- V s)f|| (12) 

and similarly for xs — xj. It can be shown that the 
exponential terms are given by 

IkR — rs + (VR — Vs)r|| + |ks— rT + fVs — V t )t|| 

-|kR-/-T + (V R -V T )r|| 

= b(xR,xs,xj) ||/> s — /> T — V ;< r || 2 /2 

where, 

v„ = (— -)(V R -V T )-(Vs-V T ) (13) 

\XR-XjJ 
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For the non-exponential terms in \jr\ (/*r(/ 2)), using Eq. (48) 
gives the same result as for \j/\ (r^(t\)), i.e., Eq. (44). 

Hence, we can write Eq. (46) at two different times as 

k 2 f 

^\(XR,p K (t\))^ — I b(x R.JC S| .*t) 
ikb(xR, xs^xj) 


exp 


/°Si ~ P T 


s'(rs ] (t\))drs ] 


(14) 


and 


k z [ 

>Ai(-*R./°r(* 2 )) % — I b(xR,xs 2 ,x T ) 
ikb{. xr,xs 2 ,-vt) 


exp 


s'(.r Sl (t2))drs 2 


Ps 2 ~ Pt ~ V I * T I 


(15) 


where the subscripts “1” and “2” on xs and p s have been 
used to indicate that we are considering scatterers at the 
points rsj and r s 2 , respectively. 

The important distinction between deterministic and ran- 
dom quantities is that the latter should be quantified via their 
statistical properties. For our purposes, we will consider 
second-order correlation functions, and subsequently the as- 
sociated Fourier frequency spectrum. Since x//\ = xi + iS l, 
it can be seen that xi = fa = (i/q -hV r *)/2 and S\ = 
Im\l/\ = (xj/] — i/c*) /2/.HBBhvhat do Im stands for? It is 
straightforward to show that the correlation function of the 
log-amplitude field is given by 


5 Xi(/°r( ? |)’/°r( , 2))= (/OR(fl))V f l (Pr('2))) 

+Re(i / 1 (/ 0 R (f|))iAf (/>r(? 2 >))] 


(16) 


For the correlation function of the phase, the first term 
in Eq. (16) changes sign; and for the cross-correlation of 
the log-amplitude and phase, the real parts are replaced by 
imaginary parts and the second term in Eq. ( 1 6) changes sign. 

In order to obtain simplified expressions for the double in- 
tegrals in Eq. (16), two further approximations are invoked. 
First, we apply Taylor’s Hypothesis (Monin and Yaglom, 
1 972) to the field s' (r s 2 fo)) = s' (r s 2 (t \ ) + Vs 2 t) . This hy- 
pothesis assumes that due to pure advection, s' (r s 2 (^)) = 
e ' (r s 2 (^i )) . That is, the permittivity fluctuation due to a scat- 
terer at rs 2 (q) will have the same numerical value after it 
moves to r s 2 (/2)- This means that the correlation function of 
the permittivity field can be written as 


B e ' (rsj (^l)i **s 2 (^2)) ={£' [rS] (/*s 2 (*2))) 
i=(s'(rs 1 (.ti))e'(rs 2 (,t\))) 


(17) 


(it is assumed that the permittivity field is a real-valued one). 
This allows for Eqs. (51) and (52) to be written solely as 
a function of t\ 9 so we can drop that notation except for 


where it is needed explicitly. Next, we invoke the Markov 
approximation (Tatarskii, 1971, Section 64). Consider the 
correlation function of the permittivity field 

5 f'(x S| -*s 2 . Ps ] . Ps 2 ) = (e'(*Si . / 0 Sl )e'(*s 2 . Ps 2 )) (! 8 ) 

For a homogeneous field, the correlation function can only be 
a function of p S] — p So and xsj and furthermore if the 

field is also isotropic it can only be a function of || p S] — p$ 2 1| 
and |^S] — xs 2 \- Under the Markov approximation it is as- 
sumed that there is no correlation of the permittivity field in 
the longitudinal direction, and hence 

^'(•*S,.JCS 2 ,/>S,./Os 2 ) = 5 (| j CSi -*%|)M|p Sl — PS*!) ( 19 ) 

where A E t is the correlation function in the trans- 
verse coordinates. Using Eq. (19), the two quantities, 

(V'iCo R ( f i))V'lO°R( f 2))) and ('Al (/M'l))^ 0 °r(' 2))) can 
be evaluated as 


(*Al (/ 0 r(M ))^1 (PR(t 2 ))) 
k 2 7T 2 


f 

i\\q II 2 ' 

/ exp 

kb(ri) 


•MMI lltWH 


( 20 ) 


<Pe' (0,||?||)|l9l|d||?||d>7 
and, 

(Vm (prM)^(pr(>2))) 

= MM\V ll (ri)[r)q>AO,M)Mdmdfl (21) 

In Eqs. (20) and (21), (p n >{ 0, ||^||) is the three-dimensional 
spectrum of the (assumed) isotropic permittivity fluctu- 
ations with q = (q y ,q z ), rj = (x s , +*s 2 )/2, and = 
( y? - JCT )^ R _ y? ) > where R=x r — xj. We need the real and imag- 
inary parts of Eqs. (20) and (21) to calculate the correlations. 
Equation (21) is real- valued, so its imaginary part is zero. 
The real and imaginary parts of Eq. (20) will result in a co- 
sine and a sine term, respectively, in the integral. Except 
for these trigonometric terms, the integrands for Eqs. (20) 
and (21) are the same. Therefore, using the trigonometric 
identity, 2sin 2 # = 1 — cos2 #, the correlation function of the 
log-amplitude fluctuations is given by 


*R oo 


&x\ (A*r(* 1 )>/ > R (^l + r )) = 4n 2 k~ j 1 9v(0.||f() 


XT 0 

\ \ ■ 2 ( Il9ll 2 ' 
t sin 


•A) (ll9ll I V^rj) ||r)sin-^ 2 ^'^ ) j \\q\\d\\q\\dr] (22) 

where we have used the spectrum of the index of refraction 
field, { p n t , instead of that for the permittivity field. (For weak 
scattering, cp F > ^4 (p n >.) We have also indicated the explicit 
ranges of the integrals. If we set Vr = Vj = 0, this means 
thatV M -> — Vs, and as above, if we move the origin to the 
transmitter we recover Eq. ( 1 9) — ( 11) in Ishimaru (1997). So 
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again, we see that our results have the correct limiting form 
as presented in the literature. Using 2cos 2 0 = 1 + cos20, 
the correlation function of the phase fluctuations is the same 
as that for the log-amplitude case, excepting that the sine- 
squared term becomes a cosine-squared term. The cross- 
correlation function is of the same form excepting that the 
sine-squared term becomes a sine term - without the factor of 
two in the denominator of the argument - and the coefficient 
in front is given by Tc 2 k 2 / 2. (After we performed the deriva- 
tion described above, a reference was found that provides a 
very similar result: Strohbehn, 1974) 

To compute the frequency spectrum, we take the cosine 
transform of Eq. (22). We define the cosine transform of a 
function g(t) as 


Amplitude Spectrum for Different 77 Values 



00 

g(f) = 4 J g(t)cos(2nft)dt (23) 

0 


Hence, 


<p X] (f) = 16n 2 k 2 J ?y(0,M)Jo(ll<7ll ||V M ||r) 

| cos(27T fr) \\q 1| d \\q \\dr)di (24) 


snr 


/ llgll 2 N 

\2kb(r}) ^ 


After performing the integration over rand ||#||we get 
(Comman, et al., 2009) 


r)]+Ari/2 


^xi(/)= f 


1 6ir 5 ' 2 k 2 AC 2 (n) 


4/3 


rj]-Ar)/2 

f2TO^)_ 

1 5r(5/6) 


VAn) 


e ia * 2 t/ [ 1 /2, — 1 /3, —iay] 


dil (25) 


where we have used the isotropic von Karman form for the 
index of refraction spectrum 


cp n '(0,\\q\\ 2 ) = AC 2 (r ] ) 


M 2 + ( 2 ”/ l) 


2~'6 


(26) 


where A is a constant, C 2 ,( 7 y) is proportional to the intensity 
of the turbulence (“structure constant”) at the point 77, R is 
the distance along the line of sight from the transmitter to the 
receiver, k is the transmitter wavenumber, and L is a length 
scale related to the correlation function for the index of re- 
fraction field. Also in Eq. (25), we are assuming a patch of 
turbulence centered at 77 = rj\ with extent A 77; and we have 
denoted a = r)(R — rj) / kR , x = 2 Ttf/V^, y = x 2 + 2tc /L 2 ; 
Tis the gamma function, and V is the confluent hypergeo- 
metric function of the second kind. The frequency spectrum 
of the phase is given by an expression similar to Eq. (25), ex- 
cepting that the minus sign in the curly brackets is replaced 
by a plus sign. If the turbulence patch is narrow enough along 
the line of sight, then the integral over 77 can be approximated 


Fig. 3. Amplitude frequency spectrum (in Hz) with varying turbu- 
lence location values. 


by a simple mid-point formula, 


<*>*>(/) 


1 6 tt 5 / 2 / c 2 A C 2 , (rj\)Ari __ 4/3 
Vpim) y 


[ 2r(i/3) 
1 5T(5/6) 


-Re [e ia * 2 U [1 /2, — 1/3, -iay] 


(27) 


where a = a(rj\ ), x =x(rj\ , /)and y = y{r ]\ ,/, L). 

Since it is a multiplicative factor, C 2 ,(rj\)Arj simply 
changes the overall amplitude of the spectrum, and changes 
in the turbulence length scale essentially raises or lowers the 
spectral amplitudes at the lowest frequencies. The changes 
in 771 are more complicated, affecting both the amplitude 
and the frequency of the oscillations - as can be seen in 
Fig. 3 between one and two Hz. The variations in these pa- 
rameters are not distinct, i.e., changes in one parameter can 
look like changes in another one. This fact makes parameter 
estimation difficult, as will be discussed subsequently. 


3 Parameter estimation 

Assuming that the received signal is a Gaussian random pro- 
cess, its spectral values at each frequency will have an expo- 
nential distribution. This result comes from the fact that the 
real and imaginary parts of the Fourier transform of the signal 
(divided by the variance) are both standard and independent 
Gaussian random variables. Hence, the spectrum, which is 
proportional to the sum of the squares of these quantities, is 
a sample from a Chi-Squared distribution with two degrees 
of freedom (cf. Theorem 8.8 in Freund, 2004). Furthermore, 
a Chi-Squared distribution with two degrees of freedom is 
an exponential distribution with mean and variance equal to 
two (cf. Definitions 6.3 and 6.4 in Freund. 2004). Therefore, 
the spectrum is distributed exponentially with mean and vari- 
ance equal to twice the variance of the signal. It is important 
to note that this randomness is due to the random permittivity 
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Amplitude Spectrum 
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Fig. 4. Amplitude frequency spectrum (in Hz) for theoretical and 
simulated case. 


field - not the additive receiver and processing noise. Using 
the fact that the spectral points are exponentially distributed, 
numerous realizations of the spectra can be generated by re- 
placing each spectral point with a sample from an exponen- 
tially distributed random number sequence. An example of 
this is shown in Fig. 4. The exponential probability density 
function is given by, 

P{x) = ^-e~ x lf > (28) 

P 

where is the mean value for the distribution. The prob- 
ability distribution function is given by P(x < b) = 1 — 
exp(— b/fi), so that P(x < ^0.63; and hence, more of- 

ten than not, samples from this distribution will be less than 
the mean, as can be seen in Fig. 4. For an exponential distri- 
bution, the standard deviation is equal to the mean. Hence, 
for the simulation, the spectral level at a given frequency 
is given by replacing that point by a sample from an expo- 
nential distribution with its mean value being the theoretical 
spectral value at that point (i.e., from Eq. 27). The simu- 
lation is intended to represent real-world situations; hence, 
it was assumed that the received signal would be available 
at 50 Hz, and approximately 10 s worth of data would go 
into a single spectrum. The parameters used in the simula- 
tion are: C~, Ar? = 4 x 10 -9 , L = 5 km, and rj \ = 20200 km. 
Furthermore, a noise floor, Q = 3.5 x 10 -5 , is added to the 
data. The magnitude of this level is derived from COSMIC 
GNSS-LEO data. 

Referring back to Eq. (27), the parameters to be estimated 
are Cj if {r]\)Ar], L, r)\ and Q : the turbulence intensity, the 
turbulence length scale, the location of the center of the tur- 
bulence patch from the transmitter, and the noise level. As 
the functional form of Eq. (27) is highly non-linear in L and 
rj i , a simultaneous solution for all four parameters is compli- 
cated. The preferred solution is via the maximum likelihood 


(ML) method (Freund, 2004). Consider a model of the form, 


P'i =c<pi(m,L)+Q (29) 

where refers to the frequency, c = C^,(r]\)Arj 9 and 
(Pi is the right-hand side of Eq. (27), divided by c. Consider 
a measured spectrum, denoted X/ . It is assumed that the X* 
are independent exponential random variables with expected 
values, (X,) = /x, . The likelihood function is then given by 


L = 


N x 
_ v — L 

^ M/ 

e 1=1 

~T v 

n ^ 

/=! 


(30) 


where N is the number of frequency values (it is assumed that 
by context, the reader will not confuse the use of L for the 
likelihood function as well as the turbulence length scale). 
The negative log-likelihood function is 


-\nL = J 2 


; = I 


'*«• . ' 
H - In /x/ 

M/ 


(31) 


The ML estimates for the parameters are found by solving 
the coupled set of non-linear equations, 3(— \nL)/dpj = 0, 
where the pj are the parameters, c, Q,L and rj\. Instead 
of dealing with this formidable problem, a simpler two- 
dimensional problem will be discussed to motivate the tech- 
niques. It will be assumed that the length scale L is known. 
Since the model, /x;, is linear in c and Q , a linear weighted 
least squares approach can be employed to estimate those two 
parameters. Once that is done, these values can then be used 
in a one-dimensional ML solution for rj \ . Some variations on 
this approach can also be used, and will be discussed below. 
Consider the least squares problem for c and Q. It will be the 
simultaneous solution to the coupled equations, 3 R/dpj = 0, 

N 

( Pj = {c, Q }), where R=^2 (M/ ~ X/) 2 u>/, and the weights 

i=l 

Wj are independent of c and Q. If the weights are chosen 
such that Wj =<p“ 1 , then this is similar to a Chi-Squared 
minimization problem. The solution for the parameters can 
then be solved via Cramer’s method for a 2 x 2 system of 
equations, giving 


c = 


N N 

E g-*E 

f,y=l 1 /=1 


N 


N 


E ^-*E*i 

ij= 1 1 i= 1 

N 

E 


N 


Q = 


(32) 






An alternative is to calculate Q directly from the high- 
frequency parts of X,. Denoting this estimate, Q, a 
one-dimensional weighted least squares solution then gives 

E (*- - Q) 

< = (33) 

X>; 

i=i 
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Log— Liklihood Function Varying rf 



Fig. 5. Negative log-likelihood function as a function of r)\ in km, 
with varying turbulence location. 


Log-Liklihood Function Varying L 



Fig. 6. Negative log-likelihood function as a function of r]\ in km, 
for varying turbulence length scale. 


A further alternative - the one that will be employed in the 
following - is to calculate c via a “pseudo” ML method. That 
is, we assume that the data V, is replaced by V, = X { — Q , 
and then solve for c via the ML formula, 3(— lnL)/c = 0. 
Letting /Z/ = ccpi 


a N 

-Y 

c!C^ 


/ = 1 


'Xi , , - ' 

— h ln/X; 
.Ml 


N r 


=-£ 


/= 1 L 


xr 

< Pi . 


+ N c 


(34) 


Equating this to zero, gives the desired ML-like estimate 
for c, 


c( )— 1 1 ~ Q 

C ^ ^ <Pi(n) 

Plugging this back into Eq. (31) gives 

N 

— InL(rj) = V[1 +c(?7)] + ^ln^,('?) 

1 = 1 


(35) 


(36) 


Where the explicit dependence on rj has been indicated (re- 
call that we are not including dependence on the length scale 
L at this point). Note that the last term is independent of the 
data and can thereby be computed a priori. Equation (36) 
is now a one-dimensional function and the solution for r)\ is 
given by minimizing it with respect to rj. For the results pre- 
sented below, we merely take the minimum value of Eq. (36) 
as a function of rj as the estimate for 77 1 . Once the estimate 
for r]\ is determined, this value is plugged back into Eq. (35) 
to estimate c. 

Different values of c will shift the negative log-likelihood 
curves up or down (for simplicity in terminology, the phrase 
“likelihood” will mean “negative log-likelihood”). Figures 5 
and 6 show the results of Eq. (36) when varying rj\, and 
L. Changes in r)\ and L change the likelihood function in 
more complicated fashions, not only varying the level of the 
curves, but also their shape. The flatness of tl)e likelihood 
function at the smaller r) poses a problem for the rj\ estima- 
tion - one would like a steep bowl leading down to the correct 


minimum. The harness means that variations in the likeli- 
hood function - due to the exponentially distributed random 
errors in the spectral values - can cause a false minimum to 
be smaller than one for the true minimum. A simple method 
was developed to deal with this problem: using a weighting 
function in the denominator of Eq. (35). This function was 
chosen to be solely dependent on rj. The only portions of 
the line of sight that are of interest are those that go through 
the atmosphere or at least portions of it where turbulence can 
be strong enough to disturb the signal significantly. For ex- 
ample, excluding outer regions of the ionosphere where the 
permittivity perturbations due to fluctuations in the electron 
density field can be small. Hence, the weighting function 
can be chosen such that it is close to one for those portions 
of interest along the line of sight, and drops smoothly to zero 
away from those regions. The effect of the weighting func- 
tion on the likelihood function is illustrated in Fig. 7. It can 
be seen that the likelihood function is unchanged essentially, 
between approximately 18 000 to 22 000 km. This is impor- 
tant, as the purpose here is not to bias the estimates, but to 
try to eliminate non-physical estimates. 

In order to examine the statistical properties of the c and ij \ 
estimators, 1 000 realizations of the simulated spectra were 
generated. Before discussing these results, however, some 
preliminary concepts are presented. 


3.1 Theoretical error analysis 


Since the random errors in ??iare not easily described theo- 
retically, these will be analyzed in the context of the simula- 
tion results. However, the error statistics of c can be studied 
analytically. Consider then the estimator for c as given by 
Eq. (36). Recall that this was referred to as a pseudo-ML es- 
timator. This is because the correct ML estimate is given by 
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Log-Likelihood Function: Weighting 



Fig. 7. The effect of the weighting function on the log-likelihood 
function as a function of 771 in km. 


- lnL =£ 


Xi 


+ ]n(c<pi + Q) 


using 

N r 

y 

f^lctpi + Q 

Then, c is found by solving 

Q d(— InL) _ f X, | <pj 

L (C<Pi + Q) 2 <Pi c< Pi + Q . 


dc 


(37) 


(38) 


which is now non-linear in c, thus requiring a numerical solu- 
tion. Therefore, we restrict our analysis to the “pseudo-ML” 
estimate. First, consider the expected value of c, 


l N 


<*/>-(£>) 


N 


~ N 2-*, 


m-Q 


N <pi(n 1) n — (pi(m) 

N 

cwvnj-ry — u 

= c 


__ j_y- c<Pi(rj l)+ Q ~ Q 


N U 


<Pi(rj 1) 


(39) 


Which shows that this is an unbiased estimator of c (where it 
has been assumed that Q is an unbiased estimator of (Land 
r] 1 is the true value for the parameter 77). Next, consider the 
variance of c - with the same assumptions made above, 


Var(c)=Var f ^ 


Xi 




(40) 


It can be shown that the variance for a linear combination of 
independent random variables, y, 9 is given by (Freund, 2004) 


Var \^2 a > yij=yafV2i(Yi ) 

Hence, with a / = \/ (Ncpi), Eq. (40) becomes 
w 1 Var (Xt) 


N [<Pi(rn)f 


(41) 


(42) 


Since it has been assumed that the Xj are independent, expo- 
nentially distributed random variables, their variance is equal 
to the square of their mean, and Eq. (42) becomes 

.. , , A {c<Pi(m)+Q ) 2 

Var c)= > 2 — (43) 

Note that if Q = 0, Var (c)=c 2 /N. 

Things are more complicated when 77 1 is not the true value. 
Denote this, rj \ . The estimated value for c (denoted c) is now 

N 

(D; l n 1 ) 

(44) 


{ c) = -^ <Pi(m) 


and hence the bias, b = {c) — (c) is given by 
N 


b = c 


] ly> P/tal) 

l) 


(45) 


Therefore, it can be seen that the bias is a function of the 
sum of the ratios of the spectrum evaluated at the true and 
estimated rj \ values, respectively. Note that if r)\ = rj\, then 
the bias will be zero, as expected. 

The variance of cis given by Eq. (43), but with fj\ in the 
denominator. Note that if Q = 0, the variance of c is 


v “®-(s ) 2 E 

i = 1 


Vital) ' 

l). 


(46) 


Since the Xj are exponentially distributed random variables, 
the random error is significant - as can be seen in Fig. 4. In 
order to reduce the random error, averaging spectra is help- 
ful. Consider the statistics of the estimator of c with aver- 
aged spectra (denoted c). The averaged spectral values are 
given by 

1 M 

*■=«£<*■>> < 47 > 


7 = 1 


where j indexes the elements that go into the average - of 
which there are M. The estimator, c, is then given by 


1 


N M ( y \ 
v A i ) j 

1 = 1 T 1 = 1 7 = 1 T 




(48) 


Interchanging the order of summation and summing over i 

M 

gives, c — Yl c j/M . Taking the expected value gives 
;=i 


1 a. 




(49) 


7=1 


where it was assumed that the (c 7 ) are identical. This shows 
that the estimator from the averaged spectra is unbiased. 
Using Eq. (41), the variance of c is given by 


Var(4-V„(|>,/*).|;^^ 


(50) 
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Distribution of C 2 n Estimates 



No Weighting Weighted No Weighting, Weighted, 


Averaged Averaged 

Fig. 8. Distribution of C^jAtj estimates over 1000 realization. The 
dashed horizontal line is the true value. 


where it was assumed that the variances of cj are identical. 
This shows that averaging the spectra reduces the variance of 
c by a factor of 1/M. 

As seen in Eqs. (44) and (46), the ratio <Pi(r)\)/<pi(rj\) is 
important in determining the error properties of the estima- 
tor c. In order to study this further, we will need the follow- 
ing concepts. The probability density function for a gamma 
distribution is given by 


P(x) = 


X ct-\ e -x/p 

P a r(a) 


(51) 


If X is gamma-distributed with parameters a and ft, and visa 
positive constant, then vX is gamma-distributed with param- 
eters a and v/3. Since the spectral values X/ are assumed to 
be independent exponentially distributed random variables - 
or gamma-distributed with a = 1 and = /x, — it follows that 
the ratio X,//z/ will be exponentially distributed with mean 
one (use the property just discussed, with v = //,). If the /x/ 
are evaluated with the estimates for Q, c and //i , then the ratio 
should be approximately distributed as an exponential distri- 
bution with mean one. The same holds for averaging the 
spectra, but here the ratio Xi/fif will be gamma-distributed 
with a = M and ft = 1/M (Freund, 2004) (recall that X, is 
given by Eq. 47). These properties can then be used to test 
the quality of the estimation methods via a goodness-of-fit 
analysis. 


Distribution of ij Estimates 



No Weighting Weighted No Weighting, Weighted, 

Averaged Averaged 


Fig. 9. Distribution of t]\ estimates. The dashed horizontal line 
indicates the true value. 


4 Simulation studies 

As mentioned above, 1000 independent realizations of the 
amplitude spectrum were generated. Equation (37) was used 
to calculate estimates of rj \ , and then those are in turn are 
used to estimate c via Eq. (36), with ri = r]\. Recall that 
the true values for these parameters are 77 1 =20 200 and 
c — 4x 10 -9 . There are four different estimation method- 
ologies that were employed: the pseudo-ML method, one in 
which a weighting function is used, and then the same meth- 
ods applied to averaged spectra. Five spectra were used in 
the averaged-spectrum analysis. A discussion of the overall 
estimator performance is presented first, followed by a set of 
case study analyzes that deal with more specific aspects. 

4.1 Overall statistical performance 

Figures 8 and 9 illustrate the overall simulation results. 
These are hybrid plots, with both the coloring and widths be- 
ing related to the number of samples at that specific param- 
eter value. Furthermore, the smaller purple boxes indicate 
outliers. Tables 1 and 2 present some descriptive statistics 
for the results. Most of the statistics used in the tables are or- 
der statistics, i.e., quantiles and inter-quantile ranges. These 
metrics are more robust in the presence of outliers, and so 
they are more appropriate here than other dispersion metrics 
such as skewness and kurtosis. It can be seen that the me- 
dian values are all quite close to the expected value, reflect- 
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ing the comments above about how the estimators are unbi- 
ased. Those comments referred to the estimator for r , but that 
fact is also represented in the median value of the estimator 
for r]\. The use of the weighting function has the desired 
effect: it eliminates the most egregious outliers. However, 
from the tables and the figures, it can be seen that the de- 
crease in outliers comes at the price of increasing the errors 
in the “middle” statistics (e.g., the 40-60 and 30-70 quan- 
tile ranges). It is clear that spectral averaging is a powerful 
method for reducing the random errors - even with averaging 
only five spectra. The beneficial effects of spectral averaging 
are far greater than using a weighting function. A definitive 
bias towards larger outliers in c can be seen in Fig. 8 . This is 
more clearly exemplified in Fig. 10, a contour-scatterplot of 
the c estimates versus the 771 estimates. The horizontal and 
vertical lines indicate the true value of those parameters. The 
solid black curve comes from the sum in Eq. (44) and subse- 
quently embodied in the equation for the bias, Eq. (45). This 
plot show how the bias is a fundamental aspect of the prob- 
lem, not an artifact of the processing methodology. Recall 
that the bias comes from the sum of the ratios of the true am- 
plitude spectrum (i.e., with the true value of ij\) divided by 
the amplitude spectrum evaluated at a different r\ \. As can 
be seen from Fig. 3, excepting for the lowest frequencies, 
the spectrum for the r 7 i = 19200 curve is always lower than 
that for r ] i = 20200 ; hence, the sum will be greater than one, 
thus producing an overestimate of c. Alternatively, if the es- 
timated value of rj\ is greater than the true value, the sum of 
the ratios will be less than one and will give an underestimate 
of c. This explains the biased structure for the estimates of c, 
as seen in Figs. 8 and 10. 

A significant concern in any estimation algorithm is the 
identification and elimination (or mitigation) of outliers. In 
the results just presented, it could be seen that the overall sta- 
tistical performance of the estimators is quite good, but in a 
real world application one wants to know if a given parame- 
ter estimate is accurate or not. This is a difficult problem in 
general and even more so in this application because the pa- 
rameters do not have independent effects on the estimators. 
As mentioned above, the ratios of the true (in this case, simu- 
lated) spectra to the model spectra with the estimated param- 
eters should be approximately distributed as an exponential 
random variable with mean one. The collection of ratios in 
question is taken over the set of frequencies samples for the 
given realization. An investigation was made into whether 
fitting an exponential distribution to the ratios could be used 
to determine the accuracy of the parameter estimates. This 
is a standard technique in the category of “goodness-of-fit” 
tests. A number of tests were employed: Anderson-Darling, 
Cramer-von Mises, Kolmogorov-Smirnov, Kuiper, Pearson 
chi-squared, and the Watson U 2 test (Lindgren, 1968). Each 
test produces a so-called test statistic, which then leads to the 
probability (or “p-value”) that the given data comes from the 
specified distribution. For each realization, the ratios were 
calculated and then tested as to how likely they were to come 


Scatter Plot of Cl Estimates vs. 7/ Estimates 



18000 19000 20000 21000 22000 

Estimated rf Values 

Fig. 10. Contour-scatterplot of C 2 Ai) estimates versus ij\ esti- 
mates. The horizontal and vertical lines indicate the true value val- 
ues for those parameters, respectively. The solid black curve is from 
Eq. (45). 

from a mean-one exponential distribution. Figure 1 1 illus- 
trates the results using the chi-squared test, giving the p-value 
as a function of the c estimate for each of the 1000 real- 
izations. The color-coding is proportional to the number of 
counts at that (c, p-value) grid cell. An optimal result would 
produce something like an upside-down “V” pattern (though 
not uniform in the number of counts), with the base of the 
two legs at the zero p-value axis, and the apex at the true c 
value and a p-value of one. That is, the higher the p-value 
(i.e., closer to one), the better the estimate - and vice-versa. 
As can be seen from the figure, the chi-squared test is clearly 
not optimal. There is skill, but there are too many mid to low 
p-values associated with good c estimates, as well as some 
poorer c estimates with high p-values. In other words, the 
p-value for this test could not be used as a single criterion 
for identifying outliers. The Pearson chi-squared test is used 
so frequently that a comment is warranted as to one reason 
that it does not perform well in this case. As discussed in 
Lindgren (1968, p. 486), the Pearson test is not well-suited 
to continuous distributions, the difficulty laying in choosing 
the appropriate binning strategy. The other tests, such as the 
Kolmogorov-Smirnov and Anderson-Darling, which use the 
distribution functions directly, do not have this limitation. 
However, it turned out that none of the individual tests had 
sufficient skill to be used in outlier detection. Thus, a sim- 
ple hybrid method was used as the goodness of fit metric - 
the maximum p-value over all the tests. The results of this 
method for the c parameter estimation are shown in Fig. 12. 
While certainly not optimal, this approach is a significant im- 
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Table 1. Statistics for maximum likelihood estimates of C~ Arj. 


C n 2 Ar? 

Median 

(xlO -9 ) 

STD 
(x 10 -9 ) 

40-60 

Quantile 

Range 

30-70 

Quantile 

Range 

20-80 

Quantile 

Range 

10-90 

Quantile 

Range 

Basic 

3.98 

0.90 

0.31 

0.67 

1.13 

2.00 

Weighted 

4.01 

070 

0.29 

0.60 

0.96 

1.64 

Spectral 

Average 

4.00 

0.24 

0.09 

0.18 

0.32 

0.50 

Weighted, 

Spectral 

Average 

4.00 

0.22 

0.08 

0.18 

0.29 

0.48 


Table 2. Statistics for maximum likelihood estimates of r /\ . 


n i 

(km) 

Median 

STD 

40-60 

Quantile 

Range 

30-70 

Quantile 

Range 

20-80 

Quantile 

Range 

10-90 

Quantile 

Range 

Basic 

20200 

903 

15 

400 

1010 

2045 

Weighted 

20 190 

653 

90 

360 

845 

1560 

Spectral 

Average 

20200 

197 

10 

40 

120 

280 

Weighted, 

Spectral 

Average 

20200 

176 

30 

60 

100 

270 


provement over the individual tests. For the r] \ estimates, 
there was more skill than in the c estimation scenario, but as 
a stand-alone outlier detection method, this goodness-of-fit 
method is still not adequate. In the next section, a number 
of case studies are discussed, which leads to a better un- 
derstanding of the difficulties in the estimation and outlier 
identification problem. 

4.2 Simulation case studies 

Five cases out of the 1000 realizations are analyzed in the 
subsequent material: realizations 64, 177, 306, 356, and 
127. The first four were chosen to study the relationship 
between p-values and parameter estimates. Specifically, the 
four cases illustrate the combinations of high and low p- 
values with both good and bad parameter estimates, (i.e., 
high-good, high-bad, low-good and low-bad). Realization 
127 is used to show how statistical comparisons between the 
log-likelihood functions derived from the simulated data and 
parameter estimates can be misleading. 

The first case studied is realization number 64. Figure 13 
through Fig. 15 show the spectra, likelihood function and dis- 


tribution of ratios for this realization, respectively. Figure 13 
shows the true spectrum (black), simulated spectrum (red), 
and the spectrum with the estimated c and rj \ values (green). 
Another spectrum function is plotted, (labeled “fit,” in blue), 
where the parameters were set “by-eye” to get a better fit to 
the likelihood function (more on this below). Figure 14 illus- 
trates the likelihood functions derived from the true spectrum 
(black); the simulated spectrum (red); the spectrum from the 
estimated parameters (green); and the spectrum from the fit- 
ted parameters (blue). Figure 15 shows a histogram of the 
spectral ratios using the spectrum derived from the estimated 
parameters, a curve representing a smoothed histogram (red), 
and a curve which is a mean-one exponential density func- 
tion (black). Table 3 gives the estimated values of c and r )\ , 
a “low” and “high” rating of the p-values, and a chi-squared 
type of error function, computed as 




|-lnL D (;?,)-(-lnL E (q,))l 2 

lnZ. E (/?,) 


(52) 


where, the sum is over all 77 values (15 000 to 25 000 km with 
a resolution of 10 km, i.e., P= 1000), -lnLo is the negative 
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Chi-Square Test 



l.xlO' 9 2.X10- 9 3.X10" 9 4xl0* 9 5.X10- 9 6.X10" 9 7.X 10- 9 8 xl0" 9 


Estimated c 2 „ Values 

Fig. 11. Probability values as a function of C^jAr) estimates: Chi- 
squared test. 


Max over all Tests 



l.xlO" 9 2.X10- 9 3.X10" 9 4.x 10- 9 5.X10- 9 6.x 10- 9 7.X10- 9 8.XIO- 9 


Estimated c\ Values 

Fig. 12. Probability values as a function of C“ Aij estimates: maxi- 
mum over all tests. 


log-likelihood calculated from the data, and -IiiLe is that 
derived from the estimated parameters. Another statistic, £f, 
is calculated using the log-likelihood function from the fitted 
parameters instead of that from the estimated ones. Note that 
these statistics are not goodness of fit metrics for the fit to the 


Amplitude Spectrum, Realization 64 



Fig. 13. Amplitude frequency spectrum (in Hz) for realization 
64, showing fit ( rj = 20 200, L = 5 km, and C^jArj = 3.84 x 10 -9 ), 
compared to estimated (rj = 17 210km, L = 5km, and C^Ar) = 
8.4 x 10" 9 ). 


Log- Likelihood Function, Realization 64 



Fig. 14. Log-likelihood functions, as a function of r]\ in km for 
realization 64, showing hand-fitted as well as estimated and true 
functions. 

spectrum, but rather for the fit to the likelihood function. The 
reason for this choice is that the parameter estimates come 
from the likelihood function - not the spectrum. 

From Table 3, it can be seen that the parameter estimates 
are quite poor for realization 64, and the error statistics are 
moderately large (relative to the other realizations). This im- 
plies that in this case, neither the Ee and E? statistics, nor the 
difference between them would indicate that the parameters 
are significantly in error. The p-values for this realization are 
very high, which indicates that they are also not helpful in 
determining whether the parameters are outliers. This is re- 
inforced by how well the distribution of the ratios matches 
the exponential curve - as seen in Fig. 15. As mentioned 
previously, the random errors for an exponential distribution 
are large, and hence the log-likelihood function can have un- 
dulations, which can result in a global minimum offset from 
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Table 3. Parameter estimates and statistical measures for selected realizations. 


Realization 

Number 

n\ 

C^rj 

P-values 

Log-Likelihood 
Error Statistic 
(£e) 

Log- Likelihood 
Error Statistic 

(£f) 

64 

17210 

8.4 x 10 -9 

High 

100 

54 

127 

20200 

4.0 x 10 -9 

High 

45 

8 

177 

20200 

4.0 x 10 -9 

Low 

16 

2 

356 

17 320 

9.1 x 10" 9 

Low 

284 

10 


Distribution Plot. Realization 64 



Fig. 15. Distributions plots of the ratios in Eq. (45), for realization 
number 64. 

the “true” minimum. It was also mentioned that this issue is 
especially problematic for smaller 77 values, where the like- 
lihood function has a shallow incline. This is exactly what 
occurs with this realization, as clearly seen in Fig. 14. The 
model parameters chosen “by eye” produce a somewhat bet- 
ter fit to the empirical log-likelihood function. This was a 
surrogate for a more sophisticated algorithm which might 
find the best fit to the empirical log-likelihood by varying 
the c and 771 (as well as L) parameters as part of a three- 
dimensional minimization algorithm. This gives hope for 
more accurate parameter estimation by using a more sophis- 
ticated minimization method, e.g., estimating L in addition 
to c and 77 1 . Since the values for c and 771 are more important 
for practical applications, the error induced by lowering the 
value for L, in exchange for more accurate c and 771 estimates 
is acceptable. However, as seen in Figs. 5, 6 , and Eq. (37), the 
effect on the log-likelihood function due to variations in these 
parameters are not independent. This is quite apparent in the 
effect that changing 77 (cf. Fig. 5) has on the log-likelihood 
function. That is, both the level and minimum value change 
with varying 77 . 

Realizations 356 and 177, both have low p-values with 
poor and excellent parameter estimates, respectively. £e is 
quite large for realization 356; however, adjusting c in a mi- 


Log- Likelihood Function, Realization 177 



Fig. 16. Log- likelihood functions, as a function of 771 in km for 
realization 177, showing hand-fitted as well as estimated and true 
functions. 


nor fashion (to 4.3 x 10 -9 ) and L in a more significant way 
(to 9), produced an excellent fit to the empirical likelihood 
function (£p = 10). Another interesting aspect of this case 
is that the fitted log-likelihood function, which produced an 
excellent fit to the empirical one, is still off in level to the true 
one. This is due to the random errors that produce a simu- 
lated log-likelihood function that is off in level. This means 
that minimizing the £p function can still produce parameters 
that are in error. Realization 177 resulted in excellent pa- 
rameter estimates, as well as a very good £e value. The fit 
did not change the parameter estimates appreciably, although 
£f is better than £e. This case is interesting in that the esti- 
mated log-likelihood function is well matched to the true one 
- as seen in Fig. 16 - whereas the fit log- likelihood function 
is well matched to the simulated one. Of course, in a real 
world scenario, the true log-likelihood function would not be 
available so one would have to depend on the fit to give all 
the information. 

This case also brings up an interesting point, one that is 
even more apparent in realization 127 (cf. Fig. 17): that good 
parameters estimates only require that the estimates of the 
level (c) and location of the minimum (771 ) match those of the 
true data, i.e., at a single point in the (rj\, — \nLv) domain. 
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Log- Likelihood Function, Realization 127 



Fig. 17. Log-likelihood functions, as a function of r]\ in km for 
realization 127, showing hand-fitted as well as estimated and true 
functions. 


On the other hand, getting a good fit to the log-likelihood 
function requires reasonable matches at many points. No 
combination of parameters might produce a good fit to the 
“measured” (real or simulated) log-likelihood function; how- 
ever, if the fit and true curve matched at just the right point, 
the estimates would be perfect. This implies that the fitting 
method should be more robust in parameter estimation - or 
at a minimum it could provide information regarding data 
quality. The hope that p-values could provide some measure 
of confidence in the parameter estimation was not fulfilled. 
Nevertheless, it was found that distributions with numerous 
outliers (i.e., large ratios) tended to produce lower p-values. 
This leads to two possible augmentations for the future: re- 
calculating the p-values with the outliers removed, and us- 
ing that as a revised confidence measure; and removing the 
outlier points from the parameter estimation process. 

4.3 Weighting and spectral averaging 

As mentioned above, other techniques that can improve pa- 
rameter estimation entails the use of a weighting function 
in calculating the log-likelihood function, and spectral av- 
eraging. Figure 18 shows the benefit in using a weighting 
function. The red and brown curves are the un- weighted 
empirical and estimated likelihood curves, respectively, for 
realization 64. Recall that the parameter estimates for this 
case were poor (c = 8.4x 10“ 9 and ij\ = 17210). Since 
the rj i estimate is on the low end of the tj— range, the use 
of a weighting function should improve the estimates. As 
can be seen from the figure (blue and green curves for the 
weighted empirical and estimated likelihood curves, respec- 
tively), the new estimates are superior to the un-weighted 
method: c = 3.98 x 10 -9 and ij\ =20150. Tables 1 and 2 
give the statistical results over 1000 realizations. Next, we 
turn to spectral averaging. A small number of realizations (5) 
are averaged to show the benefits of this method. Of course, 


Log-Liklihood Function, Weighting: Realization 64 



Fig. 18. Log-likelihood functions (as a function of r]\ in km) - 
empirical and estimated - for the basic and weighted parameter es- 
timation methods. The data are from realization 64. 

for a GNSS-LEO occultation geometry, the vertical extent of 
the atmosphere covered over the period required to average 
five spectra can be large. Running averages could be used, 
but then the samples would not be independent and the ef- 
fect of averaging on the variance - as given by Eq. (50) - 
would be lessened. Figure 19 illustrates the benefit of aver- 
aging on reducing the random errors. The red curve is the 
averaged spectrum and the black curve is the true spectrum. 
Comparing this case with, for example, the case shown in 
Fig. 4, it can be seen that the averaged spectrum is matched 
better to the true spectrum, and hence the parameter esti- 
mation should be improved similarly. Figure 20 shows the 
log-likelihood functions from the true and averaged spec- 
trum. Again, the beneficial effects of averaging can be seen, 
i.e., the log- likelihood function generated from the averaged 
spectrum is better matched to the true log-likelihood func- 
tion. Tables 1 and 2 give the statistical results over 1000 re- 
alizations. It can be seen that by reducing the random error, 
spectral averaging improves the parameter estimation much 
more significantly than the weighting scheme. It was shown 
above that when averaging spectra, the distribution function 
for the ratios becomes a gamma distribution with parameters 
ol — M and f$= \/M\ where M is the number of spectra go- 
ing into the average (here, five). This can be seen in Fig. 21 
with the simulated data. 

5 GPS-COSMIC occultation case studies 

The theory, techniques, and observations presented in the 
previous sections would be inadequate if they did not per- 
tain directly to real data. In this section, an analysis of two 
case studies from GPS-COSMIC occultations is presented. 
The purpose of these case studies is to show how the meth- 
ods and results from above apply to the real world scenario 
- not as an in-depth analysis of the real data. These cases 
were selected since they were from a region in which there 
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Distribution Plot: Spectral Averaging 
Realization 25 



Fig. 19. Averaged and theoretical frequency spectra (in Hz). Five 
spectra, beginning at realization 25, were averaged. 


Fig. 21. Distribution of ratios for the five averaged amplitude fre- 
quency spectra shown in Fig. 19, above. Note that the theoretical 
density function is from a gamma distribution. 


Log-Liklihood Function: Spectral Averaging 


Realization 25 



Fig. 20. Log-likelihood function, as a function of r]\ in km, derived 
from averaged and theoretical amplitude frequency spectra. Five 
spectra, beginning at realization 25, were averaged. 

was vigorous convective activity, as well as turbulence re- 
ports from pilots in the vicinity. In each case, a 500-sample 
analysis window was selected and processed using the same 
methodologies as used with the simulated data. There is one 
difference, however. A trend was removed from the ampli- 
tude time series data prior to the analysis. The trend was 
computed via a running median filter with a width of 2001 
points, or approximately 40 s. The appropriateness of the fil- 
ter window was checked by verifying that the residuals had a 
mean that was close to zero. The noise level was computed 
by taking the median value of the log-amplitude spectrum 
over the last 50 points (i.e., the high-frequency region). 

Figure 22 shows the amplitude time series over the course 
of occultation 2 (arbitrary nomenclature). The black curve is 
the raw data, the red curve shows the trend, and the two ver- 


SNR Data - Trend Removal 
COSMIC Occultation 2, (2300,2800) 



Fig. 22. SNR data for COSMIC occultation number 2 (black). 
Trend curve is in red and the vertical, dashed blue lines indicate the 
analysis window. 


tical dashed blue lines delineate the analysis window. This 
window was chosen because it had large-amplitude fluctua- 
tions, and could be used to test subjectively whether some of 
the assumptions presented in Sect. 2.1 were violated; specifi- 
cally whether inhomogeneity and/or layered structures would 
affect the results. Figure 23 shows the spectra of the data 
(red), one using the estimated parameters (in green, with 
c = 5. 1 x 1CT 10 , rn = 19010, L = 5.0 and Q= 7.0 x 10“ 6 ), 
and one using the fitted parameters (in blue, with c = 5. 19 x 
10~'°, tj\ = 19010, L = 3.1 and Q = 6.0 x 10~ 6 ). The like- 
lihood functions for the data (red), the estimated parameters 
(green) and fitted parameters (blue) are presented in Fig. 24. 
For this case, £ E = 31.6 and £f = 2.1. Figure 25 illus- 
trates the distribution of ratios. From these three figures it 
can be seen that the characteristics of the real data are sim- 
ilar to those of the simulated cases (e.g., realization 64). 
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Fig. 23. Amplitude frequency spectrum (in Hz) for occultation 2 
(red), model spectrum using estimated parameters (green), and “by- 
eye” fitted parameters (blue). 


Fig. 24. Log-likelihood functions, as a function of ij\ in km, 
from data (red), estimated parameters (green), and fitted parame- 
ters (blue) for occultation 1 . 

The p-values for this case are mid-level to high, as might 
be expected by how well the distribution from the real data 
matches an exponential one. 

The amplitude time series for occultation 1 is shown in 
Fig. 26. It can be seen that there are significant fluctuations 
within the window. This window was chosen because of that 
fact - to see how well the theory applies to such a dynamic 
case. Figure 27 shows the spectra for the data (red), the spec- 
tra derived from the estimated parameters (green), and that 
using a set of fitted parameters (blue). The spectrum of the 
data shows features that were not seen in the simulated data. 
For example, the “hump” from around 3 to 6Flz. In fact, 
if we look at the amplitude (red) and phase (black) time se- 
ries over the analysis window (Fig. 28), a significant ramp 


Distribution Plot 


Cosmic Occultation 2, (2300,2800) 



Fig. 25. Distribution of ratios for occultation 2. 


can be seen in the amplitude data at around 54s, with what 
appears to be a damped oscillation in the subsequent 1.5 s. 
The phase data shows a step-like structure during the same 
1.5 s, with the steps in concert with those in the amplitude 
data. (It should be noted that the de-trending used with the 
phase data was very different than that used with the am- 
plitude data. A polynomial fit over an extended time seg- 
ment, centered at the analysis window was employed. For 
this example, the fit was performed over a window extend- 
ing 150 points before and after the analysis window.) The 
source of these oscillations is unknown, but their periods are 
approximately 10-16 samples, which at 50 Hz. gives fre- 
quencies from around 3-5 Hz - which matches what is seen 
in the frequency spectrum. The hump in the spectrum gen- 
erates magnified oscillations in the likelihood function, seen 
in the empirical (red) curve in Fig. 29, with the unfortunate 
by-product that one of these oscillations happens to be the 
global minimum, in turn skewing the parameter estimation 
(see the green curve in this figure). A by-eye fit was per- 
formed (blue curve) which provided far superior parameter 
estimates - or at least a better fit to the empirical likelihood 
function - since the true values are unknown. The estimated 
and fitted error statistics are Ee = 299.2 and Ee = 2.4, re- 
spectively. For the estimated parameters, c = 2.76x 10 -9 , 
rj\ = 17730, Q = 1 .8 x 10 -5 and L = 5, and for the fitted pa- 
rameters, c= 1.31 x 10 -9 , ij\ =20800, Q =6.0x 10~ 6 ,and 
L =7.0. It is encouraging that even with the data quality 
issues, physically reasonable parameter estimates can be ob- 
tained with fits to the likelihood function. The distribution of 
the ratios for this data also contains artifacts from the oscil- 
lations. In fact, these are directly related to the hump in the 
spectrum around 5 Hz, as those values - being much larger 
than the model spectrum - in turn produce large ratios. These 
large ratios in turn produce poor p-values. 
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SNR Data - Trend Remov al 


COSMIC Occultation 1 , (2450,2950) 



Fig. 26. SNR data for COSMIC occultation number 1 (black). 
Trend curve is in red and the vertical, dashed blue lines indicate 
the analysis window. 


Fig. 27. Amplitude frequency spectrum (in Hz) for occultation I 
(red), model spectrum using estimated parameters (green), and “by- 
eye” fitted parameters (blue). 

6 Conclusions 

In this paper, an overview of the derivation used to produce 
the model spectrum was presented. This derivation started 
from the Helmholtz wave equation, and using a reasonable 
set of assumptions, resulted in the model spectrum, Eq. (26), 
and its mid-point approximation, Eq. (28). This latter model 
was used in the subsequent analysis. The assumptions used 
in the derivation were spelled out clearly. It was pointed out 
that, as we considered this a first-order solution, two real 
world atmospheric situations, inhomogeneity (including lay- 
ered structures) and anisotropy, were not considered. The 
model derivation was based on standard weak scattering the- 
ory which can be found in many references. However, in our 
case we needed to modify the standard methods to include 
the motions of the transmitter, receiver, and atmosphere, as 
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Amplitude and Phase Data (Trend Removed) 
COSMIC Occultation 1.(2450,2950) 



time (s) 

Fig. 28. De-trended amplitude and phase time series (in seconds) 
from the analysis window for occultation 1 . (Between the vertical, 
dashed blue lines in Fig. 26.) Note spike just after 54s, and oscilla- 
tions (“ringing”) for the subsequent 1 .5 s. 


Fig. 29. Log-likelihood functions, as a function of rj j in km, 
from data (red), estimated parameters (green), and fitted parame- 
ters (blue) for occultation 2. 

well as placing the origin of the coordinate system at an arbi- 
trary point. It was shown that our results match those in the 
literature when making the appropriate simplifications. 

A parameter estimation methodology was presented, fol- 
lowed by a discussion of the statistical properties of the esti- 
mators. A detailed study using simulated data was presented 
after the theoretical discussions. A number of metrics were 
examined to evaluate the accuracy of the estimates, e.g., me- 
dian, standard deviation, and various inter-quantile ranges. It 
was found that the estimates of the turbulence intensity and 
the location of the turbulence along the line of sight could 
be retrieved in an accurate manner. In this instance, the term 
“accurate” refers to the global statistics over the 1000 simu- 
lated realizations. In a practical application, one would like 
to have accurate estimates for each measurement, or at least 
a notion as to when the estimates are unreliable. Goodness 
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of fit metrics (p-values) were used to see if the distribution 
of the ratios of the measured/simulated spectrum values to 
those from the true one could be used as an indicator of poor 
estimates. It was found that this technique was not reliable 
in that regard. Four cases were shown that demonstrated the 
problem, one where the p-values were high and the estimates 
were poor, and vice versa - as well as cases where the p- 
value was high and the estimates were good, and vice versa. 
This was also borne out with a scatterplot of the parameter 
estimates and p-values. Even when choosing the maximum 
p-value over the different goodness of fit tests, the correla- 
tion of p-value and quality of the parameter estimate was 
reasonable - but not great. 

Since the parameter estimation technique used the global 
minimum of the log-likelihood function, random errors in the 
spectrum could result in significant errors when trying to es- 
timate the location of the turbulence along the line of sight. 
An important step in improving the parameter estimation was 
fitting the model-based likelihood function to the empiri- 
cal one. In addition to incorporating the turbulence length 
scale parameter into the fitting, a weighting function was in- 
troduced into the likelihood estimation methodology. The 
weights were close to one for most of the ij values deemed 
to be in the earth’s atmosphere, and decreased rapidly out- 
side those bounds. This not only allowed for the elimination 
of the most egregious outliers, but it also lead to improved 
parameter estimates from the data interior to the weighting 
function’s effective region. Finally, it was shown that, due 
to reducing the random errors, spectral averaging greatly im- 
proved parameter estimation. These smaller errors then prop- 
agate through the parameter estimation resulting in more ac- 
curate estimates. It was seem that spectral averaging was far 
superior to the weighting function in improving the parame- 
ter estimates. The methods developed for the analysis of the 
simulated data was then applied to the analysis of real occul- 
tation data. In the two cases presented, it could be seen that 
the general characteristics seen in the simulated data were 
also seen in the real data - which is encouraging. 

Future work includes the implementation of a three- 
parameter estimation algorithm including the length scale, 
and a more thorough analysis of real occultation data. This 
is especially needed for phase data. Some analysis has been 
performed with simulated and real phase data, with encour- 
aging results, but the extensive analysis presented here with 
the amplitude data has not been performed. There are dis- 
tinct issues regarding trend removal with the phase time 
series. A local polynomial fitting method was discussed 
briefly; however, the frequency content in the data resulting 
from polynomial fitting is highly dependent on the order of 
the polynomial used and the window over which the fit is 
performed. This indicates that a more robust methodology 
needs to be developed. 

As mentioned above, future work should include the ef- 
fects of inhomogeneities (including layered structures) and 
anisotropy. The reason that the inclusion of layered media is 


difficult is because large-scale layers will refract the incom- 
ing wave, which then means that the straight-line approxima- 
tion for the propagation path needs to be modified. The per- 
mittivity field is separated into background and random parts. 
(Note that in this context, “background” refers to all aspects 
of the deterministic permittivity field.) One could then re- 
write the original differential equation in a non-Cartesian co- 
ordinate system, e.g., a local one which has as its axes the 
tangent, normal and binormal vectors to the “main” propa- 
gation path - for example that described by geometrical op- 
tics. These are so-called trajectory coordinates (see for exam- 
ple: Hill, 1985; Mazur and Felsen, 1987). Another technique 
is to use a path integral approach to determine the Green’s 
function in a multiple scattering context (Mazur, 2002). Yet 
another approach is to use the so-called Distorted- Wave Born 
(or Rytov) Approximations (see for example: Beylkin and 
Oristaglio, 1985; Devaney, 1979). In this method, one as- 
sumes that the permittivity field is separated into background 
and random parts, however the heuristic model is that the in- 
coming waves are refracted by the background and are then 
scattered by the random field. The contribution of all such 
scattered fields is then what is measured at the receiver. That 
is, it is a single-scattering process, but the electric field at the 
scatterer is now “distorted” by the background field. Of the 
three methods, this is probably the most tractable one. 
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